rm(list = ls())
set.seed(123)

# Load Packages

library(tidyverse) # 1.3.1
library(stargazer) # 5.2.2
library(scales) # 1.1.1
library(modelsummary) # 0.7.0
library(vtable) # 1.3.1
library(gtsummary) # 1.4.1
library(gridExtra) # 2.3
library(jtools) # 2.1.3
library(sandwich) # 3.0-1
library(lmtest) # 0.9-38
library(sensemakr) # 0.1.3
library(MatchIt) # 4.2.0

# load data
data <- read.csv("perspective_taking_merged.csv")

## Filter out those who refused informed consent 
# + failed attention checks

data %>% 
  filter(informed_consent == 1 & attention_check1 == 3) %>%
  mutate_all(na_if,"") %>%
  filter(age >= 18) %>%
  filter(!is.na(attention_check1)) -> data

#
SEs <- function(data){
  empty_data <- c()
  data1 <- vcovHC(data, type = "HC1")
  empty_data <- sqrt(diag(data1))
  return(empty_data)
}

#### Add Lived Experience indicator
# Plot COVID experience
data %>%
  mutate(covid_ever = ifelse(covid_contracted == 1 | ## Label individuals with lived COVID-19 experience 
                               covid_self == 1, "Yes", "No")) %>%
  mutate(treatment2 = factor(treatment, ## Arrange the treatment in the right order
                             levels = c("Control Undoc",
                                        "Treat Undoc",
                                        "Cancer",
                                        "Immigrant") ),
         pid_3level = recode(pid_7level, # Make the 3 level party ID
                             `1` = 1, ## Democrats
                             `2` = 1,
                             `3` = 1,
                             `4` = 2, ## Independents
                             `5` = 3, ## Republicans
                             `6` = 3,
                             `7` = 3),
         polid_3level = recode(pol_ideology, # Make the 3 level pol id
                               `1` = 1, ## Liberals
                               `2` = 1,
                               `3` = 1,
                               `4` = 2, ## Moderates
                               `5` = 3, ## Conservatives
                               `6` = 3,
                               `7` = 3)) -> data


########################################
############ Table A15 ##############
########################################     

# Focus on control and COVID-19 treatment 
livedexp <- data %>% 
  filter(treatment %in% c("Control Undoc", 
                          "Treat Undoc") ) %>%
  droplevels()

# Run unadjusted regression
livedexp_unadj <- lm(response_deserving ~ treatment*covid_ever, livedexp)
livedexp_unadj_se <- as.data.frame(SEs(livedexp_unadj))

# Run adjusted regression
livedexp_adj <- lm(response_deserving ~ treatment*covid_ever
                     + age + gender + ethnicity + as.factor(pid_3level) + educ + income, livedexp)
livedexp_adj_se <- as.data.frame(SEs(livedexp_adj))

## create modelsummary
modelsummary(list(
  "Lived Experience\nUnadjusted" = livedexp_unadj,
  "Lived Experience\nAdjusted" = livedexp_adj),
  stars = T, 
  vcov = "HC1",
  coef_rename = c("treatmentTreat Undoc" = "Perspective Taking\nUndocumented",
                  "covid_everYes" = "COVID-19 Experience",
                  "treatmentTreat Undoc × covid_everYes" = "Treatment COVID-19 \n + Undocumented \n x COVID19 Experience"),
  out = "default")

## Create stargazer
stargazer(livedexp_unadj,
          livedexp_adj,
          type = "text", align = T,
          font.size = "footnotesize",
          digits = 3, 
          no.space = T,
          initial.zero = F,
          keep = c("Constant", "treatmentTreat Undoc",
                   "covid_everYes", "treatmentTreat Undoc × covid_everYes"),
          se = list(livedexp_unadj_se[,1],
                    livedexp_adj_se[,1]),
          column.labels = c("Unadjusted", "Adjusted"),
          model.numbers = F,
          covariate.labels = c("Undocumented + Treatment COVID-19",
                               "COVID-19 Lived Experience",
                               "Undocumented + Treatment COVID-19 x COVID19 Lived Experience",
                               "Constant"),
          dep.var.labels = "Undocumented Immigrants as Deserving of Gov Healthcare")



########################################
############## Figure A11 ##############
########################################     


### Plot the results
plot_summs(list("Unadjusted" = lm(response_deserving ~ treatment2*covid_ever, data),
                "Adjusted" = lm(response_deserving ~ treatment2*covid_ever
                                + age + gender + ethnicity + as.factor(pid_3level) + educ + income, data) ),
           coefs = c("Treatment Cancer \n + Undocumented" = "treatment2Cancer",
                     "Treatment COVID-19 \n + Undocumented" = "treatment2Treat Undoc",
                     "Treatment COVID-19 \n + Immigrant" = "treatment2Immigrant",
                     "COVID-19 Experience" = "covid_everYes",
                     "Treatment Cancer \n + Undocumented \n + COVID19 Experience" = "treatment2Cancer:covid_everYes",
                     "Treatment COVID-19 \n + Undocumented \n + COVID19 Experience" = "treatment2Treat Undoc:covid_everYes",
                     "Treatment COVID-19 \n + Immigrant \n + COVID19 Experience" = "treatment2Immigrant:covid_everYes"),
           robust = T, inner_ci_level = 0.9,
           colors = c("dark blue", "light blue")) +
  theme(legend.position = "top") 
ggsave("appendix_figure_a11_final.jpeg")

### Tabular results
modelsummary(list("Unadjusted" = lm(response_deserving ~ treatment2*covid_ever, data),
                  "Adjusted" = lm(response_deserving ~ treatment2*covid_ever
                                  + age + gender + ethnicity + as.factor(pid_3level) + educ + income, data) ),
             coefs = c("Treatment Cancer \n + Undocumented" = "treatment2Cancer",
                       "Treatment COVID-19 \n + Undocumented" = "treatment2Treat Undoc",
                       "Treatment COVID-19 \n + Immigrant" = "treatment2Immigrant",
                       "COVID-19 Experience" = "covid_everYes",
                       "Treatment Cancer \n + Undocumented \n + COVID19 Experience" = "treatment2Cancer:covid_everYes",
                       "Treatment COVID-19 \n + Undocumented \n + COVID19 Experience" = "treatment2Treat Undoc:covid_everYes",
                       "Treatment COVID-19 \n + Immigrant \n + COVID19 Experience" = "treatment2Immigrant:covid_everYes"),
             robust = T,
             stars = T,
             vcov = "HC1",
             coef_rename = c("treatment2Treat Undoc"= "Treatment COVID-19 + Undocumented",  
                             "treatment2Cancer" = "Treatment Cancer + Undocumented", 
                             "treatment2Immigrant" = "Treatment COVID-19 + Immigrant",
                             "covid_everYes" = "COVID-19 Experience",
                             "treatment2Treat Undoc × covid_everYes" = "Treatment COVID-19",
                             "treatment2Cancer:covid_everYes" = "Treatment Cancer \n + Undocumented \n + COVID19 Experience",
                             "treatment2Treat Undoc:covid_everYes" = "Treatment COVID-19 \n + Undocumented \n + COVID19 Experience",
                             "treatment2Immigrant:covid_everYes" = "Treatment COVID-19 \n + Immigrant \n + COVID19 Experience"),
             out = "default")




###########################################
##### PROPENSITY SCORE  ###################
###########################################

# Prepare data 
data %>%
  filter(!is.na(covid_ever), # Need to remove the NAs on key variables - COVID-19 ever
         !is.na(age), # Remove NAs in Age
         !is.na(educ), # Remove NAs in Education
         !is.na(income), # Remove NAs in Income
         !is.na(treatment2), # Remove NAs in treatment
         !is.na(ethnicity), # Remove NAs in Ethnicity
         !is.na(pid_3level)) %>% # Remove NAs in Party ID
  filter(treatment2 %in% c("Control Undoc", # Remove Low Income groups
                          "Treat Undoc",
                          "Cancer",
                          "Immigrant")) %>%
  mutate(covid_ever = recode(covid_ever, # Label the COVID-19 Lived Exp
                             "No" = 0,
                             "Yes" = 1),
         treatment2 = factor(treatment2, # Arrange in the right order
                             levels = c("Control Undoc",
                                        "Treat Undoc",
                                        "Cancer",
                                        "Immigrant"))) -> pscore

### Check balances prior to matching
# Get variable list
vars_pscore <- c('age', 'educ', 'income', 'ethnicity', 'pid_3level')

# 
pscore %>%
  group_by(covid_ever) %>%
  select(one_of(vars_pscore)) %>%
  summarise_all(funs(mean(., na.rm = T)))

##
m.out0 <- matchit(covid_ever ~ age 
                  + income + ethnicity 
                  + educ + pid_3level, 
                  data = pscore,
                  method = NULL, distance = "glm")

summary(m.out0)
plot(summary(m.out0) ) ## if it does not work - input dev.off() then try again


# matching
m.out1 <- matchit(covid_ever ~ age 
                  + income + ethnicity
                  + educ + pid_3level, 
                  data = pscore,
                  method = "nearest", distance = "glm",
                  ratio = 1)
summary(m.out1) ## if it does not work - input dev.off() then try again

########################################
############## Figure A12 ##############
########################################     

# If the figure does not show up
# input dev.off() then try again

plot(summary(m.out1))
ggsave("appendix_figure_a12_final.jpeg")

#########################
#########################

## Testing regression
m.data1 <- match.data(m.out1)
m.data0 <- match.data(m.out0)


# Check balance pre and post propensity score
# First shows post-score, second shows pre-score
# Post
m.data1 %>%
  group_by(covid_ever) %>%
  select(one_of(vars_pscore)) %>%
  summarise_all(funs(mean)) -> post; post

# Pre
m.data0 %>%
  group_by(covid_ever) %>%
  select(one_of(vars_pscore)) %>%
  summarise_all(funs(mean)) -> pre; pre


### Run t tests to make sure they are balanced

# Party ID
t.test(m.data1$pid_3level[m.data1$covid_ever == 1],
       m.data1$pid_3level[m.data1$covid_ever == 0],
       var.equal = F)

# Age
t.test(m.data1$age[m.data1$covid_ever == 1],
       m.data1$age[m.data1$covid_ever == 0],
       var.equal = F)

# Education
t.test(m.data1$educ[m.data1$covid_ever == 1],
       m.data1$educ[m.data1$covid_ever == 0],
       var.equal = F)

# Income
t.test(m.data1$income[m.data1$covid_ever == 1],
       m.data1$income[m.data1$covid_ever == 0],
       var.equal = F)

# Ethnicity
t.test(m.data1$ethnicity[m.data1$covid_ever == 1],
       m.data1$ethnicity[m.data1$covid_ever == 0],
       var.equal = F)

## Testing regression
fit1 <- lm(response_deserving ~ treatment2*covid_ever 
           + age + income + ethnicity + educ
           + pid_3level, 
           data = m.data1, weights = weights) 

coeftest(fit1, vcov. = vcovCL, cluster = ~subclass) 



######## Pre vs post score matching plot

# remove low inc groups
pscore %>%
  filter(!treatment %in% c("Control Low Inc",
                           "Treat Low Inc")) -> pscore

# Make a vector of coefficients
variables = c("treatment2Treat Undoc" = "Treatment Undocumented \n COVID-19",
              "treatment2Cancer" = "Treatment Undocumented \n  Cancer",
              "treatment2Immigrant" = "Treatment Immigrant \n COVID-19",
              "covid_ever" = "COVID-19 Yes",
              "age" = "Age",
              "educ" = "Education",
              "income" = "Income",
              "ethnicity" = "Ethnicity",
              "pid_3level" = "PID 3 Level",
              "treatment2Treat Undoc:covid_ever" = "Treatment Undocumented \n COVID-19 x COVID-19 Yes",
              "treatment2Cancer:covid_ever" = "Treatment Undocumented \n Cancer x COVID-19 Yes",
              "treatment2Immigrant:covid_ever" = "Treatment \n Immigrant x COVID-19 Yes")

########################################
############## Table A16 ##############
########################################     

modelsummary(list("Not Matched" = lm(response_deserving ~ treatment2*covid_ever
                                     + age + educ + income + ethnicity 
                                     + as.factor(pid_3level), 
                                     data = pscore),
                  "Matched" = lm(response_deserving ~ treatment2*covid_ever
                                 + age + educ + income + ethnicity 
                                 + as.factor(pid_3level), 
                                 data = m.data1, weights = weights)),
             vcov = "robust", 
             stars = T,
             coef_omit = "Intercept",
             coef_map = variables)


########################################
######## Sensitivity Tests #############
########################################     

########################################
####### Table A17 and Figure A13 #######
########################################     

data %>%
  filter(treatment2 %in% c("Control Undoc", # Keep only treatment with undocumented 
                           "Treat Undoc")) %>%
  mutate(treatment2 = case_when(treatment2 == "Control Undoc" ~ 0, # label treatment
                                treatment2 == "Treat Undoc" ~ 1),
         covid_ever = case_when(covid_ever == "No" ~ 0, # Label COVID-19 lived experience
                                covid_ever == "Yes" ~ 1)) -> study_sens

## Adjusted model
model.sens <- lm(response_deserving ~ covid_ever
                 + age + educ + income + gender + ethnicity + pid_3level, study_sens)

# Use sensemakr
covid.sensitivity <- sensemakr(model = model.sens, 
                               treatment = "covid_ever",
                               benchmark_covariates = "pid_3level",
                               kd = c(1:3) )

########################################
############# Table A16  ###############
########################################

### Table results for COVID Ever
ovb_minimal_reporting(covid.sensitivity, format = "latex")

########################################
########### Figure A13  ###############
########################################

plot(covid.sensitivity)
ggsave("appendix_figure_a13_final.jpeg")





########################################
############# Figure A14 ###############
########################################     

# If the figure does not show up
# use dev.off() then try again

#### COVID-19 Severity
data %>%
  mutate(treatment2 = factor(treatment, # Arrange the treatment in a better order
                             levels = c("Control Undoc",
                                        "Treat Undoc",
                                        "Cancer",
                                        "Immigrant",
                                        "Control Low Inc",
                                        "Treat Low Inc"))) %>%
  group_by(covid_symptoms, treatment2) %>%
  summarise(sample_count = n(),
            mean = mean(response_deserving, na.rm = T),
            se = sd(response_deserving) / sqrt(n())) %>%
  filter(!is.na(covid_symptoms),
         !is.na(treatment2)) %>%
  ggplot(aes(x = covid_symptoms, y = mean,
             ymin = mean - (2*se),
             ymax = mean + (2*se)) ) +
  geom_pointrange(position = position_dodge2(w = 0.5)) +
  labs(x = "") +
  facet_wrap(treatment2 ~.,
             nrow = 3,
             labeller = as_labeller(c(`Cancer` = "Treatment \nUndocumented \nCancer",
                                      `Control Undoc` = "Control \nUndocumented \nCOVID-19",
                                      `Treat Undoc` = "Treatment \nUndocumented \nCOVID-19",
                                      `Immigrant` = "Treatment \nImmigrant \nCOVID-19",
                                      `Control Low Inc` = "Control \nLow Income \nCOVID-19",
                                      `Treat Low Inc` = "Treatment \nLow Income \nCOVID-19")) ) +
  scale_y_continuous(labels = scales::percent) +
  labs(x = "COVID-19 Severity", title = "", y = "Mean Deservingness Support") +
  scale_x_continuous(breaks = seq(0,5, by = 1)) +
  theme_bw() + 
  theme(legend.position = "top") 
ggsave("appendix_figure_a14_final.jpeg")



########################################
############# Figure A15 ###############
########################################     
data %>%
  mutate(treatment2 = factor(treatment,
                             levels = c("Control Undoc",
                                        "Treat Undoc",
                                        "Cancer",
                                        "Immigrant",
                                        "Control Low Inc",
                                        "Treat Low Inc"))) %>%
  filter(treatment2 %in% c("Control Undoc", "Treat Undoc")) %>%
  group_by(covid_symptoms, treatment2) %>%
  summarise(sample_count = n(),
            mean = mean(response_deserving, na.rm = T),
            se = sd(response_deserving) / sqrt(n())) %>%
  filter(!is.na(covid_symptoms),
         !is.na(treatment2)) %>%
  ggplot(aes(x = covid_symptoms, y = mean,
             ymin = mean - (2*se),
             ymax = mean + (2*se),
             shape = as.factor(treatment2),
             color = as.factor(treatment2)) ) +
  geom_pointrange(position = position_dodge2(w = 0.5)) +
  labs(x = "") + 
  scale_y_continuous(labels = scales::percent) +
  scale_color_manual(name = "",
                     labels = c("Control",
                                "Perspective Taking"),
                     values = c("blue",
                                "red3")) +
  scale_shape_discrete(name = "",
                       labels = c("Control",
                                  "Perspective Taking")) +
  labs(x = "COVID-19 Severity", title = "", y = "Mean Deservingness Support") +
  scale_x_continuous(breaks = seq(0,5, by = 1)) +
  theme_bw() + 
  theme(legend.position = "top", 
        axis.text.x = element_text(size = rel(2)),
        axis.text.y = element_text(size = rel(1.5)),
        axis.title.x = element_text(size = rel(1.55)),
        strip.text = element_text(size = 15),
        legend.key.size = unit(1, "cm"),
        legend.title = element_text(size = 11),
        legend.text = element_text(size=15)) 
ggsave("appendix_figure_a15_final.jpeg")




########################################
############# Table A18 ################
########################################

data %>%
  mutate(ethnicity = factor(ethnicity),
         income = factor(case_when(
           income == 1 ~ "20k or less",
           income >= 2 & income <= 6 ~ "20 - 59k",
           income >= 7 & income <= 9 ~ "60 - 99k",
           income == 12 ~ "100k or more",
           TRUE ~ "N/A"),
           levels = c("20k or less", "20 - 59k",
                      "60 - 99k", "100k or more",
                      "N/A")),
         region = factor(case_when(
           region == 1 ~ "Northeast",
           region == 2 ~ "Midwest",
           region == 3 ~ "South",
           region == 4 ~ "West")),
         gender = factor(case_when(gender == 1 ~ "Male",
                            gender == 2 ~ "Women")),
         educ = factor(case_when(
           educ == 1 ~ "Less than HS", 
           educ == 2 ~ "HS",
           educ == 3 ~ "Some college",
           educ == 4 ~ "AA",
           educ == 5 ~ "BA",
           educ == 6 ~ "Postgraduate",
           TRUE ~ "N/A"),
           levels = c("Less than HS", "HS",
                      "Some college", "AA",
                      "BA", "Postgraduate", "N/A")),
         age = factor(case_when(
           age >= 18 & age <= 29 ~ "18 - 29",
           age >= 30 & age <= 41 ~ "30 - 41",
           age >= 42 & age <= 54 ~ "42 - 54",
           age >= 55 ~ "55 or more",
           TRUE ~ "N/A")),
         ethnicity = factor(case_when(
           ethnicity == 1 ~ "White",
           ethnicity == 2 ~ "Hispanic",
           ethnicity == 3 ~ "Black",
           ethnicity == 4 ~ "Asian",
           ethnicity == 5 ~ "Native American",
           TRUE ~ "other")),
         pid_3level = factor(case_when(
           pid_3level == 1 ~ "Democrat",
           pid_3level == 2 ~ "Independent",
           pid_3level == 3 ~ "Republican")),
         polid_3level = factor(case_when(
           polid_3level == 1 ~ "Liberal",
           polid_3level == 2 ~ "Moderate",
           polid_3level == 3 ~ "Conservative")) ) %>%
  select(covid_symptoms, gender, educ, age,
         region, ethnicity, income,
         pid_3level, polid_3level) %>%
  sumtable(group = "covid_symptoms")


